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Abstract 

We present a detailed study of a model for strain-induced metal-insulator phase coexistence in perovskite 
manganites. Both nanoscale and mesoscale inhomogeneities are self-consistently described using atomic 
scale modes and their associated constraint equations. We also examine the stability of domain config- 
urations against uniform and nonuniform modifications of domain walls. Our results show that the long 
range interactions between strain fields and the complex energy landscape with multiple metastable states 
play essential roles in stabilizing metal-insulator phase coexistence, as observed in perovskite manganites. 
We elaborate on the modes, constraint equations, energies, and energy gradients that form the basis of our 
simulation results. 
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I. INTRODUCTION 



Over the last several years, much attention has focused on the multiscale inhomogeneities ob- 
served in perovskite manganites ^ Unlike inhomogeneities for some other complex electronic sys- 
tems, such as stripes in high-T c cuprates, the coexistence of metallic and insulating phases within 
the same crystals of manganites has been directly observed through various high resolution local 
probes, such as dark field images and scanning microscopy. 2 4 Nanoscale inhomogeneities have 
also been implied based on x-ray diffraction results. 5 Theories based on chemical randomness and 
electronic phase separation have been proposed as a mechanism for such inhomogeneities. 6 7 How- 
ever, theories based on chemical randomness alone assume an exact degeneracy between metallic 
and insulating phases, 7 which ultimately leads to a homogeneous phased Also, the effect of the 
Coulomb interaction has not been incorporated adequately into the model describing electronic 
phase separation. 6 

We have previously proposed an intrinsic mechanism for phase coexistence,^ in which the 
interaction between strain fields plays an important role, as speculated in earlier literature 
Specifically, our model 9 includes intrinsic complexity of the energy landscape and long range 
anisotropic interaction between strain fields, and shows how such physics can give rise to multi- 
scale inhomogeneities observed in manganites. Our theoretical idea is supported by experimental 
results.^ For example, the observed large scale inhomogeneity of the order of 10 fim without any 
observable chemical inhomogeneity at a length scale of 0.5 fim or larger, suggests an intrinsic 
mechanism for the phase coexistence. 14 The lamellar morphology of coexisting phases observed 
in manganites and the change of domain configurations upon thermal cycling between 10 K and 
300 K further support this point of view. 15 It is also found that the anisotropic epitaxial strain in 
thin films gives rise to anisotropic percolation, which suggests that the origin of phase coexistence 
is much more strongly affected by long range strain rather than by local chemical inhomogeneity 
due to doping.EO 

In contrast to the point of view in Refs. [6] and [7] that considers extrinsic mechanism such as 
chemical inhomogeneity and disorder as being key to understanding coexistence of metallic and 
insulating phases, our work proposes that the intrinsic mechanism, that is, the competition between 
the short and long range interactions, creates a delicate energy landscape that leads to domain for- 
mation. Further, the domain walls can be pinned by the atomistic Peierls-Nabarro force, rather than 
by disorder which is frequently attributed as being the cause for the phenomenon. 6 7 We proposed 
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the basic ideas in Ref. HI but now develop them further and describe more extensive implemen- 
tation of the approach here. Specifically, we present details of our model, methods, and results, 
as well as further simulations on the stability of phase coexistence. In Sect. II, the details of the 
Hamiltonian used for the simulations in Ref. 9, and results obtained with various initial conditions 
and parameters, are presented. We also contrast these results with simulations for a system that 
include short range interactions only. In Sect. Ill, we discuss the mechanism underlying the sta- 
bility of micron-scale phase coexistence through further simulations and analysis. In Sect. IV, we 
summarize our main results. Expressions for energies and energy gradients used for simulations 
of the inhomogeneous states are provided in Appendices A and B, respectively. 

II. MODEL FOR STRAIN-INDUCED METAL-INSULATOR PHASE COEXISTENCE 
A. Properties of manganites and requirements for phase coexistence in manganites 

Perovskite manganites typically have the chemical formula REi_ x AK x Wh\Oz, where RE rep- 
resents rare earth elements, such as La, Nd, and Pr, and AK represents alkali metal elements such 
as Ca and SrJ 1 ' 17 * 18 * The important electrons for both electronic properties and structures are the 
e g electrons on Mn ions: the degeneracy of the e g orbitals leads to a strong Jahn-Teller electron- 
lattice coupling. Shortly after the discovery of colossal magnetoresistance in these materials,^ 
the importance of strong electron-lattice coupling was pointed out to explain the large resistivity 
above the ferromagnetic transition temperature in terms of dynamic Jahn-Teller polarons. 18 The 
same strong electron-lattice coupling is also responsible for a large structural difference between 
the low temperature metallic and insulating phases of these materials. In the insulating phase, the 
e g electrons are localized and the charge density forms an ordered pattern. The orbital states of the 
e g electrons, which are linear combinations of two e g orbital states, x 2 — y 2 and 3z 2 — r 2 , are also 
ordered. Due to the static Jahn-Teller effect, such a charge and orbital ordered state accompanies 
uniform and short wavelength lattice distortions. For example, the long Mn-0 bond (or the elon- 
gated e g electron orbital) in Lao.sCao.sMnOs alternates its direction in a plane, which gives rise 
to short wavelength lattice distortions. 20 Along the direction perpendicular to this plane, the short 
Mn-0 bond repeats itself, and therefore the unit cell is compressed along this orientation, giving 
uniform or long- wavelength lattice distortions. In contrast, the lattice in the metallic ground state 
of manganites has a structure close to an ideal cubic perovskite structure without Jahn-Teller lattice 
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distortions because the electrons are delocalized. The ground state can be changed between metal- 
lic and insulating phases in various ways, such as by the size of RE and AK, applied magnetic 
fields, or applied pressurespHSH 

We propose that the first key to the understanding of the phase coexistence in manganites 
is the metastability, which has been observed in many experiments. For example, the distorted 
insulating phase of manganites can be transformed into the undistorted metallic phase by either 
x-rays^ or magnetic fields,^ and the metallic phase does not revert to the insulating phase even 
after the external perturbation is removed. In particular, in x-ray experiments, 24 the reduction of 
the superlattice peak intensity and the simultaneous increase of conductivity, while the sample 
is exposed to the x-rays, demonstrate the transformation of the insulating phase into the metallic 
phase and the presence of inhomogeneity. However, an energy landscape with local and global 
energy minima is not sufficient to explain the observed sub-micron scale inhomogeneity, because 
such inhomogeneity is unstable against thermal fluctuations. As pointed out in Ref.fTTl an unusual 
aspect of inhomogeneity in manganites is its stability over a 100 K range in temperature, which is 
an indication of an extra mechanism at play that affects phase coexistence in manganites. 

Based on the strong electron-lattice coupling mentioned above and experiments showing the 
important role of strain in metal-insulator transition in manganites, 25 we propose that the extra 
mechanism should be related to the long range anisotropic interaction between strain fields. It is 
thus essential to consider the energy landscape in terms of the lattice distortion variables, which 
will ultimately have a bearing on other degrees of freedom such as magnetic moment or elec- 
tron density. The origin of the long range anisotropic interaction within this framework is the 
bonding constraint, often referred to as strain compatibility. The compatibility condition enforces 
single-valued strain fields without broken bondsP^The anisotropy reflects the discrete rotational 
symmetry associated with the lattice structures. Such long-range anisotropic interactions are re- 
sponsible, for example, for well-defined structural twin boundaries^ over distances of 100 \xm. 

The cubic phase of perovskite manganites contains five atoms per unit cell. The insulating 
charge and orbital ordered phase consists of a zig-zag pattern of the long Mn-0 bond orientation, 
which further increases the number of atoms per unit cell. Inclusion of such details is necessary for 
a complete description of properties of these materials. For the current study, however, we wish 
to focus on the following three key features of manganites essential for multiscale inhomogeneity, 
and capture them in a simple model. First, the metallic phase has almost no lattice distortions in 
comparison with the charge and orbital ordered insulating phase. Second, the insulating ground 
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state has a uniform or long wavelength (k ~ 0) lattice distortion. This property is essential be- 
cause it is the long wavelength distortion, not the short wavelength one, that gives rise to the 
long range anisotropic interaction between strain fields. Third, the insulating phase has a short 
wavelength lattice distortion, in addition to the uniform distortion. As we will show below, the 
symmetry-allowed coupling between uniform and short- wavelength distortions gives rise naturally 
to an energy landscape with multiple minima. 

B. Model system, variables and constraint equations 

Before we introduce our model, we examine whether a simpler 2-dimensional (2D) model can 
be used instead of a 3-dimensional (3D) model, in particular to capture the effect of the long- 
range strain- strain interactions. In D-dimensional space (D = 2 or 3), the anisotropic strain- strain 
interaction decays as l/r D , where r represents the distance between two points P^ED The S p a . 
tial integration of l/r D would give rise to a logarithmic divergence in both 2 and 3 dimensional 
space, 31 which indicates that the effect of the interaction would be similar for both cases. Indeed, 
recent simulations of strains in 2 and 3 dimensional space show very similar resultsP^Thus, we 
limit ourselves in this work to a 2D model for simplicity. 

One of the simplest lattices in 2D space is the square lattice with a monatomic basis shown in 
Fig.[T] By considering one isotropic electron orbital per site and nearest neighbor electron hopping, 
the lattice supports a metallic electron density of states (DOS) without a gap. Therefore, such an 
undistorted square lattice shares the first property for manganites mentioned above. To include the 
second property, we deform the square unit cell of the lattice to a rectangular unit cell, either along 
the horizontal or vertical directions. To include the third property, we incorporate the (it, n) type 
displacements of atoms along the horizontal (vertical) direction for a rectangular lattice elongated 
along the vertical (horizontal) direction. Rectangular lattices with such short wavelength lattice 
distortions support an electron DOS with a gap at the center, if we consider the natural electron 
hopping amplitude modulation by the changes in interatomic distances. Therefore, such lattices 
have an insulating DOS for the electron density of half an electron per site. Even though the 2D 
lattice described above is simple, it shares all the three properties of manganites that we believe are 
essential for the observed inhomogeneity, and provides a testing ground for whether the complex 
energy landscape and long range strain-strain interaction can indeed give rise to self-consistent 
multiscale inhomogeneities. 
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We require an energy expression for which the undistorted and distorted lattices described 
above are the local and global energy minimum states. For this purpose, we use an atomic-scale 
mode-based description of lattice distortions that we developed recently.^ In this method, we 
use normal modes of a square plaquette of four atoms, instead of displacement variables, to de- 
scribe lattice distortions. These atomic scale modes for the monatomic square lattice are shown 
in Fig. |2j The first three modes are long wave length modes, since they can be obtained by uni- 
formly deforming the square lattice. The last two modes, which correspond to (ir, it) staggered 
distortions of the lattice, are short wavelength modes. For a square lattice, each atom is shared by 
four neighboring plaquettes, which makes the modes at neighboring plaquettes dependent on each 
other. Such a constraint can be expressed in terms of equations in the Fourier transformed space, 
and the five modes can describe any lattice distortion for the square lattice with a monatomic ba- 
sis. In the long wavelength limit, the three long wavelength modes become identical to the familiar 
strain modes, which makes our approach ideal for describing nano-and micro-meter scale inhomo- 
geneities within the same theoretical framework. The inclusion of constraints allows our method 
to automatically generate the effects of the long range anisotropic interaction, the origin of which 
is the short-range bonding constraint. 

We consider an iV x N square lattice with a modified periodic condition explained below. The 
displacement variables for the atom at the site i are it? and ut. The distance between the nearest 
neighbor atoms, a, is irrelevant to our formalism presented below, and can be chosen depending 
on the relative size of the distortions compared to the lattice constants. In all the figures in this 
paper, a is chosen as 10 so that the size of the distortion relative to the interatomic distance is of the 
same order of the magnitude as observed in charge and orbital ordered manganites. In general, the 
displacement of atoms in a periodic structure can be described using two components. One is the 
component that changes monotonically as the site indices shift along a direction. This component, 
represented by a superscript 'nF' below, can not be Fourier-transformed and corresponds to the 
uniform distortion of the lattice. The rest of the displacement, represented by a superscript 'F' 
below, can be Fourier-transformed and is subject to the periodic boundary condition. Therefore, 
we express the displacements as follows: 



x x,nF i x,F 

= Mr + , 

i % i 1 

v V,nF i y.F 

U% = UV + U% , 



(1) 

(2) 
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where 



x,nF xx ■ I xy • 

y,nF xy ■ , i/y • 

W^. — fcg 2 X + fcg ly, (ft-) 



and 



U- 
i 



y,F 
UV 



J2 u l eik4 > (5) 

k 

J2 u l ei ^- (6) 



We note that u% and u y - are obtained through the Fourier transformation of u*' F and v% ,F , rather 

k k i i 

than u% and uH. The periodic boundary condition results in 

K = ^, (7) 

where = -N/2 - 1, JV/2 and n y = -N/2 - 1, ...,JV/2. For u?' nF and uf nF , the rigid 
rotation of the whole system is excluded, since it is irrelevant to the potential energy change. 
Similarly, the k = components of u^' F and u^' F correspond to the rigid translation of the whole 
system, which are set to zero. 

For the square lattice, we define the symmetry modes shown in Fig.|2]as follows: 

ei& = A (-uf -4 + «S in - ul in - ui + ul (9) 



2^2 



2V2 



— ui 

1 


+ «f+10 


— vJi 
i+10 




i+01 


^ "i+11 


— ui 

1 


- u* 

"i+10 


+ ul 
i+10 


+ u Uoi 


u 

— 

i+01 


+ u ln 


+ ul 


+ U f+10 


i+10 


"i+01 


u 

— u" 

i+01 


+ u !+ll 



i+11 



= A? ( -ti? - - + ul ln + M ? - < m + «Ln + «L, ) ■ < .10) 



i+11 



eJT) = -i- (-u% + ul + u2 1n + mH - w2 n1 -«H m +«? J . l1 -^ 1l )> (ll) 

2\/2 V * i H-10 i+10 »+01 4 +01 2+11 i+11/ 

s x (i) = 2 ( M f - wf+io " M f+oi + "f+n) > ( 12 ) 
= 2 ( M ? " M lio - M loi + M f+n) • < 13) 

These modes are fully subject to the periodic boundary condition, e.g., e\(i x , i y ) = ei(i x +N,i y ) = 
e-i(i x ,iy + N), unlike the displacement variables. Thus, they can be Fourier-transformed, for 
example, according to 

ei(i) = $>i® e<£r - (14) 
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From the definitions, we find 



ei(k = 0) 



-xx i jyy 



V2 



ex, 



-.''.V 



e 2 (k = 0) = = e 2 , 



e 3 (fc = 0) 



v/2 



= e 3 , 



S a (fc = 0) =0, 

Sy (A? = 0) = 0, 



(15) 
(16) 

(17) 

(18) 
(19) 



and 



ei(k ^ 0) 

e 2 (k^0) 

e 3 (k ^ 0) 

s x (k^0) 
s y (k^0) 



1 



-(1 - e ife *)(l + e l **)«£ - (1 + e <fc »)(l - e**»)« 



2^ L 
1 



2V2 
1 



-(1 + e ife *)(l - e lfc »)u5 - (1 - e ^)(l + e^)w| 



(1 - e ifc *)(l + e ifc »)u| + (1 + e^)(l - e ik *)ul 



2^2 

i(l-e^)(l-e^)«| 



(l-e i *«)(l- e < *»W' 



(20) 
(21) 
(22) 

(23) 
(24) 



2 V /v ' * 

We note that the k = components of the symmetry modes are from u^' nF and u y ' nF , while = 
components of u*' F and do not contribute to the distortion modes. 

The five variables are related by three constraint equations, because only two physically in- 
dependent displacement variables exist for each site. As discussed in Ref. 28, these constraint 
equations are found from the relations between the symmetry modes and the displacement vari- 
ables in reciprocal space. For k x 7^ and k y ^ 0, we invert the linear relations between [s x (k), 



s y (k)] and ui] in Eqs. (23) and (24) and replace them in the expressions with other modes in 



Eqs. (f20])-((22]). This leads to 



Js* ^ J^> J^> t Js> } 

sin — cos —s x (k) + cos — sin —sJk) — V2i sin — sin —eAk) = 0, 
2 2 v ' 2 2 yy 1 2 2 y ' 

Js> J^> # Js> J^> ^ Js> Js> j 

cos — sin —s T (k) + sin — cos —s v (k) — y2i sin — sin —eo{k) = 0, 
2 2 y 1 2 2 vK ' 2 2 y ' 

J^> ^ J^> J/i t jj-? j 

sin — cos —sJk) — cos — sin —sJk) — y2i sin — sin —e^(k) = 0. 
2 2 w 2 2 yy 1 * 2 2 w 

These constraint equations indicate that ei(7r, tt), 62(71", 7r), and 63(7?, 7r) vanish and s a 



(25) 
(26) 
(27) 

S x (7T,7r) 



and s y = s y (n, tc) are independent variables. Constraint equations for k x = or k y = should be 



considered separately from Eqs. (25l-(27|. Equations ( 15 >-( 19 > show that ei(k = 0), e 2 (k = 0) 
and e 3 (k = 0) are independent of each other, and s x (k = 0) and s y (k = 0) vanish. For k x = 



and k y ^ 0, Eqs. (20)-(24) show that e\{k) = —e 3 (k) and e 2 (k) are independent variables, and 



s x (k) = s y (k) = 0. Similarly, for k x 7^ and k y = 0, e\{k) = e 3 {k) and e 2 {k) are independent 
variables, and s x (k) = s y (k) = 0. 

To describe lattice distortions in our simulations, we primarily use the variables s x (i) and s y (i). 
These variables can be assigned arbitrarily except that they should satisfy s x (k) = s y {k) = if 



or k y = 0, as required by Eqs. (18), (19), (23), and (24). In our numerical simulations, 



we implement this condition by subtracting unphysical components with k x = or k y = from 
s x (i) and s y (i), each time we initialize or change s x (J) and s y (J). However, s x (i) and s y (J) do not 
uniquely determine lattice distortions, because of the singular relation between [s x (k), s y (k)] and 



[tt|, u|] in Eqs. (|23|) and (|24|). As seen above, ei(k = 0), e 2 (k = 0), e 3 (k = 0), ei(k x = 0, k y ^ 



0) = -e 3 {k x = 0, k y ^ 0), e 2 {k x = 0, k y ^ 0), ei(A^ 7^ 0, k y = 0) = e 3 {k x ^ 0,k y = 0), and 
e2(^x ^ Q,k y = 0) should be specified, in addition to s x (i) and s y (i), for the complete description 
of lattice distortions. 

From these variables, displacement variables, and u%, are calculated. For the non-periodic 
parts of displacements, u^' nF and ul ,nF in Eqs. (J3J) and (4), we use Eqs. ( 15 )-( 17 ) to obtain 

, nF ei(fc = 0) + e 3 (k = 0) 



ia, + A/2e 2 (A; = 0); 



ul- F = V2e 2 (k = 0)4 + ei(fc = 0) ' e3(fc = 0) z 

y/2 



v 



(28) 
(29) 



We find the periodic parts of the displacement, uv and , through the Fourier transformation 



of and u%, which are obtained by inverting two non-singular equations among Eqs. (20)-(24). 



Therefore, if k x 7^ and k y 7^ 0, we invert Eqs. (23 ) and (24 ) to obtain 



u 



k x ^0,k y ^0 



2 



U 



-S„(fc). 



(30) 
(31) 



If fca 7^ and k y = 0, Eqs. (20) and (21 ) lead to 



72 



u 



y 

k x ^0,k y =0 



1-e 
v/2 



; k -^ik). 

tn,x 



(32) 
(33) 



Similarly, if k x = and k y ^ 0, we obtain 



V2 



<=o, fc ^o = ir^eaW. ( 34 ) 

<=o, fc ^o = Y3^ ei ^)' (35) 

The k x = and k y = components of the displacements correspond to rigid displacements, which 
are set to zero: 

u% x=0jky=0 = 0, (36) 
<=o*=o = °- (37) 

By adding periodic and non-periodic parts of displacements according to Eqs. ([T|) and ([2]), we find 

u% and u%. 



C. Total energy of the model and the Hamiltonian for electronic property calculations 

In terms of the above modes, we consider the following energy expression, E tot , as the total 
energy of a model system for strain-induced phase coexistence^ 



Etot = E s + Ei + E c , (38) 

I i i , . fro ^ ^ J~f 1 „ „ f/i n „ „ „ „ 

,(39) 



A\ 2,^2 2,^3 2 

2 ei + 2 e 2 + 2 e 3 



(40) 



£ c = ^[C3(4-4)e 3 ] r ( 41 ) 

i 

The first term with short wavelength modes includes all symmetry-allowed terms up to sixth 
order with all coefficients positive, since we are interested in a first-order-transition-like energy 
landscape. The second term Ei with long wavelength modes up to second order mediates the long 
range anisotropic interactions. The third term E c represents the coupling between the long and 
the short wavelength modes, where the mode is coupled to the s x and s y modes in a symmetry- 
allowed form. This last term gives rise to the global energy minimum state with long and short 
wavelength distortions, in addition to the local energy minimum state without distortion. The en- 
ergy expression E tot gives rise to the desired energy landscape in appropriate ranges of parameters. 
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To establish a connection between the electronic properties and the lattice distortions as ob- 
served in manganites, namely metallic and insulating states for the undistorted and distorted 
phases, respectively, we use the following Su-Shrieffer-Heeger (SSH) Hamiltonian for the elec- 
tronic structure calculations in our model: 



SSH 



1 - a(v%, nn s - ui 

v i+(10) t 



4 C ?+(10) + C 



t (A 

i+(io) v 



-to l-a(t^ +(01) -i*? 



4 c ?+(oi) + 4+(oi) c ?; • 



(42) 



Here, we consider only one orbital per site and neglect the electron spin for simplicity. The oper- 
ator ct is the creation operator of an electron at a site i. In this Hamiltonian, the electron hopping 
amplitude is assumed to be linearly modified by the change in the nearest neighbor interatomic 
distances. We make the adiabatic assumption that the total energy E tot is obtained by minimiz- 
ing the energy of the system with respect to all degrees of freedom, including the electronic one, 
except for the lattice degrees of freedom. Therefore, E to t is used for the calculation of the en- 
ergy landscape and the Euler simulations, whereas H SSH is used for the calculations of electronic 
properties associated with templates of lattice distortions. 

The SSH Hamiltonian is a Hamiltonian for independent electrons, and can be diagonalized 
within a one-electron basis. Therefore, we construct electronic Hamiltonian matrices for given 
lattice distortions, u% and u%, with the basis set {cl|0)}, where |0) represents the state without 
electrons. We diagonalize the matrices numerically, and fill the eigenstates with electrons accord- 
ing to the electron density. Representing the Z-th lowest energy eigenstate as 

IO = Evtl°>' (43) 

% 

the local DOS at site % is calculated by 

D I (E) = J2^(E-E l )\z l /, (44) 
i 

which reveals local electronic properties. The same approach has been used in Ref. [32] to study 
electronic inhomogeneities around structural twin and antiphase boundaries in model systems with 
a strong electron-lattice coupling. 
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D. Energy landscape for homogeneous states 



We expect the ground state of the model energy expression E tot to be homogeneous with e\,e-i, 



§3, s x and s y distortions only, defined in Eqs. ( 15 »-( 17 ) and below Eq. (27), considering the way 



that the energy terms are selected. Therefore, we study the following energy expression, which 
includes these particular distortions only, to understand the energy landscape for the homogeneous 
states: 

Et t = E h s + E\ + E h c , (45) 

^ = f K + *1) + ^ K + Sj) + ^ + f K + Sj) + f ^1 (*1 + Sj) , (46) 
^4i ~o ^2 ~o ^4 3 ~o 

j£ = ye? + f <g + f 4 (47) 

^ = ^3 K - 3) e ~3- (48) 

Because ei, £2, S3, S x and S y are independent of each other, we minimize E^ ot with respect to ei, 
e 2 , and e 3 independently and obtain e.\ = 0, e 2 = 0, and e 3 = — C 3 (s^ — Sy)/A 3 . We insert these 
back into E^ ot and obtain the following energy expression in terms of s x and s y only: 

rph,min R -, / s~i2 \ i / fi2\ 

111 tot n f~2 , ~2\ , 1 o°3 \ /~4 , ~4\ , 1 / , o°3 \ ~2~2 



jv, 2 « + »1) + i (g. - aj) W + « + \ (g. + 2f ) #J 

+^(»1 + »1) + ^N»1(sS+»1)- < 4< » 



We find parameter values, for which E^] nm /N 2 has one local energy minimum state without dis- 
tortion and four symmetry-related degenerate global energy minimum states with distortions. Nec- 
essary conditions for such first-order-transition-like energy landscape are G[ = G± — 2C 2 /A 3 < 
and (G'1) 2 > ABH\, for which the energy minima occur at (s x , s y ) = (0, 0), (±s , 0), and (0, ±s ) 
with 



The locations of the five energy minima in the s x — s y plane are indicated in Fig.[3ja). Correspond- 
ing uniform mode distortions are e 3 = — C 3 sl/A 3 for (s x , s y ) = (±s , 0), and e 3 = C 3 sl/A 3 for 
(s x , s y ) = (0, ±sq). For comparison, we choose two sets of parameter values, one giving a shallow 
and the other a deep local energy minimum at (s x , s y ) = (0, 0), as shown with a thin blue and a 
thick red curve, respectively, in Fig.^b). In manganites, the difference in the depth of the energy 
landscape can be related to the size of rare earth or alkali metal elements, which is known experi- 
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mentally to influence the physical properties of manganites.^ Alternatively, we may consider this 
a measure of "micro strain" pi 



E. Methods of simulations for inhomogeneous states 

The energy landscape is much more complicated for inhomogeneous states because of the con- 
straint relations among the distortion modes. To study inhomogeneous configurations, particularly, 
metastable configurations, we first minimize E tot analytically with respect to all the independent 
variables except s x (i) and s y (i), that is, e\(k = 0), e 2 (k = 0), e 3 (k = 0), e\(k x = 0,^/0) = 
-e 3 (k x = 0,k y ^ 0), e 2 (k x = 0, k y ^ 0), e x {k x ^ 0, k y = 0) = e 3 {k x ^ 0, k y = 0), and 
&2{k x 7^ 0,k y = 0), and obtain an energy expression E' tot (s x , s y ). The details of the derivation and 
expression for E' tot (s x , s y ) are provided in Appendix A. 

In our simulations, we set initial configurations of s x (J) and s y (i) and relax the lattice according 
to the Euler method, 



OS x (l) 

V V dS y {i) 



(51) 
(52) 



where the superscript n or n + 1 represents the number of Euler steps taken from the initial con- 
figuration, and 7 controls the size of the Euler step. Expressions for dE' tot (s x , s y ) /ds x (i) and 
dE' tot (s x , Sy)/ds y (i) are provided in Appendix B. We change s x (i) and s y (i) for all z's simulta- 
neously at each step. We run the simulation until E' tot (s x , s y ) does not decrease further, but only 
fluctuates, which is an indication that the system has reached a local energy minimum configura- 
tion. 



F. Initial conditions and results of the simulations for inhomogeneous states 

We describe initial conditions, parameters, and results of the simulations in this subsection. 
Figure [4] shows the results of the simulations carried out on a 32 x 32 lattice for the energy land- 
scape with a shallow local minimum, shown in thin blue curve in Fig.[3jb) for homogeneous states. 
The color of each plaquette represents p 3 (i) = s x (i) 2 — s y (i) 2 , and the vertices and distortions of 
the plaquettes represent the actual locations of atoms and actual distortions. Through the coupling 
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between p 3 and e 3 in E c , positive and negative values of p 3 are usually accompanied by an e 3 
distortion elongated along y and x direction, respectively. Most plaquettes with p 3 close to zero 
have little distortion. Starting from an initial configuration of s x (i) and s y (J), randomly chosen 
between — 2s and 2s , as shown in Fig.[4ja), the system is relaxed through the Euler method with 
7 = 0.0015. Figures |4jb)-|4jh) correspond to the configurations at the Euler step n =100, 400, 
1000, 2000, 4000, and 6000, and the stable configuration at n =100000, respectivelyP 

Following the energy gradient from a random initial configuration, the simulation approxi- 
mately represents a rapid quenching of the system from a very high temperature to K. The result 
shows that most of the system initially changes into the undistorted state, as suggested in Fig.[4|cl). 
Random fluctuations of p 3 tend to cancel each other when averaged over several interatomic dis- 
tances, which prevents most regions from evolving into a global energy minimum state with lattice 
distortions. Some regions with relatively large distortions in the initial configuration evolve into 
a distorted state, as shown in Figs. |4|b)-|4jd), and nucleate the distorted phase. These distorted 
regions expand into the undistorted region [Figs.[4je)-|4jg)], and eventually transforms most of the 
system into a distorted global energy minimum state separated by anti-phase boundaries, as shown 
in Fig. |4j n )- We also find that there is a critical size for this nucleation: a distorted region of a 
small size in Fig. [4](c) changes into the undistorted state in Fig.|4jd). Such nucleation and growth 
observed in our simulations of the rapid quenching, and the presence of the critical size for the 
nucleation, reflect the first-order-transition-like energy landscape, and are features observed even 
in systems with only short range interactions. 

For contrast, we study a system with a similar energy landscape, but without the long range 
interaction. For this we consider the following energy expression: 

E m = £ yM« 2 + —M(i) 2 + -M(z) 6 

t 

Z 



2 



{[M(i+ (1, 0)) - M(z)] 2 + [M(i+ (0, 1)) - M(?)] 2 } , (53) 



where M{i) is defined at each site. The last term becomes Z(VM) 2 /2, the familiar Ginzburg- 
Landau gradient term, in the continuum limit. The gradient term of Em with respect to M{i) for 
the Euler method is 

WM{T) + IM(i) 3 + YM{i) 5 



dM(T) 

+ Z[AM(j) - M(T- (1, 0)) - M(i + (1, 0)) - M{%- (0, 1)) - M{i+ (0, 1))].(54) 
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The parameters W, X, and Y are chosen to be identical to the parameters B, G[, and Hi for the 
lattice model with a shallow local minimum in Fig. [3jb). We choose Z = 5, similar to A lf A 2 
and A 3 , since e%(i), e 2 (i), and e 3 (i) are related to the gradient of s x (i) and s y (i) in the continuum 
limit. 28 The uniform ground state for Em is M(i) = ±M , where 



-X + VX 2 - 4WY 
\/ • (55) 

The selected parameter values result in M = 0.38, identical to so for the lattice distortion model. 
Figure [5] shows the results of the simulations on a 32 x 32 lattice, in which p = M 2 , analogous to p% 
for the lattice distortion model, is shown. The initial M{i), shown in Fig. [5^ a), is chosen randomly 
between — 5.3M and 5.3M , and is relaxed according to the Euler method with 7 = 0.00015. The 
configurations at the Euler step n = 1000, 4000, 10000, 2000, 40000, and 60000 are shown in 
Figs.[5];b)-[5jg). The final stable configuration at n = 1000000 is shown in Fig.[5pi). These Euler 
steps are chosen so that they are consistent with the Euler steps in Fig.|4]after being multiplied by 
7. The system with a short range anisotropic interaction in Fig.|5]also shows nucleation and growth 
of the low energy phase. However, comparing Figs. [4] and [5]reveals distinct features present only 
in the lattice distortion model. 

First, the nucleation droplets in Figs. |4jc)-|4|e) are highly anisotropic, in contrast with those 
in Figs. [5|c)-|5je). Second, distortions separated by relatively large distance along the diagonal 
direction interact with each other, grow toward each other, and merge through the long range in- 
teraction, as seen for the yellow and red band along the 135 degree orientations in Figs.|4jc)-|4jf). 
Third, the nucleation occurs via pairs of distortions with different orientations to minimize the in- 
terface energy between the distorted region and the undistorted background, as shown in Figs.[4|c) 
and[4jd). Such features are absent in Fig. [5], where the interaction is purely short-ranged. Recent 
x-ray scattering experiments have revealed the presence of short-range anisotropic precursor cor- 
relations in the orthorhombic phase of manganites at high temperatures, which disappear in the 
rhombohedral phase. 5 Such a feature has a similarity with the anisotropic droplets observed in our 
simulations, and is reminiscent of the precursor embryonic fluctuation in martensitic transforma- 
tions. 35 The quasi-elastic central peak observed in manganite^^ near the metal-insulator phase 
transition temperature is also likely to have a structural origin, similar to the central peak observed 
in ferroelectrics. Such experimental observations 38 and similar features seen in our simulations 
indicate that the strain plays an important role in the formation of nanometer scale inhomogeneity 
in manganites. 
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To demonstrate electronic inhomogeneity associated with the structural inhomogeneity, we cal- 
culate electronic properties for the template of the lattice distortions in Fig.[4]T). We use the SSH 



Hamiltonian in Eq. (42) with t = 1 and a = 1. The typical local electron densities of states 
within undistorted and distorted regions are shown in Fig. (6fb). The local DOS is symmetric 
about E — 0, and a gap (or a "pseudogap") opens near E = in distorted regions. The small 
DOS within the gap for the distorted region is due to electron wavefunctions exponentially de- 
caying from the undistorted region. Therefore, distorted and undistorted regions have insulating 
and metallic electron DOS at a half filling without any spatial charge inhomogeneity. The map 
of the local electron DOS at E = is shown in Fig. [6^a). The possible inhomogeneity in local 
DOS without any charge inhomogeneity in our model is in contrast with other explanations for 
the inhomogeneity based on electronic phase separation, an idea similar to the phase separation in 
binary alloys. 6 

Figure [7] shows results of a simulation with parameters identical to those in Fig. [4] except for 
a narrower range of the random initial values of s x (J) and s y (i) between — s and so. Instead of 
multiple nucleations, only one nucleation emerges within the 32 x 32 lattice, which grows and 
evolves the whole system into a periodic patten of stripes with positive and negative p^. The final 
state shown in Fig.[7]T) is another metastable phase not considered in Fig. [3], where distortions only 
with a wave vector (n, tt) are considered. The result shows that multiple inequivalent metastable 
phases exist even in this simple model, and the coexistence of more than two phases is possible, 
as suggested in some manganites. 39 We use an even narrower range of the random initial s x (i) and 
s y (i) between — s /2 and s /2, in which case the system fails to nucleate a distorted region and 
remains in the undistorted phase, showing the characteristic metastability of systems with a first- 
order-transition energy landscape. The above results indicate that the low temperature metastable 
configurations depend sensitively on how the configurations are obtained, consistent with path- 
dependent experiments in manganites, such as the sensitivity to the cooling rate or strain glass 
behaviors.^ 

For the deep local energy minimum case described by the thick red curve in Fig.^b), the sim- 
ulation of rapid quenching using the Euler method for a 64 x 64 lattice does not create nucleation 
of the low energy phase. Instead, we always obtain the undistorted homogeneous state as the fi- 
nal state, which is an indication of strong metastability due to a higher energy barrier between 
the distorted and the undistorted states. In crystals, we expect line or planar defects, as well as 
thermal fluctuations, would assist nucleation. Simulations of such processes require more compu- 
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tational resources. Therefore, we start from a predesigned initial condition and relax the lattice to 
obtain stable coexistence of distorted and undistorted domains. The initial condition is chosen on 
a 64 x 64 lattice according to 



sJi x ,i v ) = s (-l) i " +H ' <^ cos 



2n(i x + i y -A) 
N 



+ 0.5}, (56) 



s y {i x ,iy) = 0, (57) 

where N = 64. The initial configuration is relaxed with 7 = 0.0002, and the stable configuration is 
obtained, which is shown in Fig.[8|a). We find stable coexistence of large undistorted [green region 
in Fig.[8|a)] and distorted [red region in Fig.[8fa)] domains, unlike the shallow local minimum case 
studied above. The size of the domain is determined only by the initial condition, and therefore 
can be as large as several micrometers, consistent with experiments for manganites. 

For comparison, we carry out simulations for a N x N system with a short range interaction 



M(i x , iy) = M < cos 



only, described by E M in Eq. (53 ) with the parameter W = 2.0, with a similar predesigned initial 
condition, 

'^^J+c,}. (58) 

For c s = 0.5, 7 = 0.0002, and iV = 64, we find that the final stable configuration is a uniform 
ground state for this system with a short interaction only, rather than a state with domains. For 
c s = 0, 7 = 0.0002, and iV = 128, we find only a line of atoms, rather than a domain, with M 
close to zero between regions with M = M and M = — M . This comparison shows that the 
strain-strain long range interaction indeed plays an essential role for the coexistence of distorted 
and undistorted phases. 

For the configuration in Fig. [8|a), the local DOS versus energy is calculated at the centers of 
the undistorted region and the distorted region, which is shown in Fig. [9} The local DOS within 
the undistorted regions shows a metallic DOS without a gap, whereas the local DOS within the 
distorted region shows an insulating DOS with a gap around E = 0. Therefore, for the chemical 
potential chosen at E = inside this gap, we obtain the coexistence of metallic undistorted and 
insulating distorted regions, as shown in Fig.[8]^b), similar to experiment results. We also find that 
the interface between the metallic and insulating regions is rather sharp, consistent with STM im- 
ages of atomically sharp interface between metallic and insulating domains in manganites. 3 Our 
results indicate that chemical inhomogeneity is not a necessary condition to have a large scale 
coexistence of metallic and insulating domains, which is in contrast to other theories. 7 Although 
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the lattice defects or segregation of dopants could play a role in nucleation, the stability of coexis- 
tence relies on the intrinsic energy landscape, which explains why external perturbations such as 
focused x-rays,^ lightj^ or electron beams alter metallic and insulating domains. 



III. STABILITY OF PHASE COEXISTENCE 



A. Stability against uniform domain wall motions 



In this section, we examine the stability of the phase coexistence against various kinds of per- 
turbations. First, we examine the energy barrier blocking a uniform shift of the domain boundaries, 
which would convert the undistorted high energy phase into the distorted low energy phase. Red 



dots connected by the lowest lines in Fig. 10 show s x (i x ,i y ) x (— l) l *+ l w versus i x for i y — 1 
near the boundary between the undistorted (i.e., i x < 51) and distorted (i.e., i x > 52) phases for 
the configuration in Fig. ^ a). To find the energy barrier against uniform domain wall shift, we 
increase the value of s x (i x , i y ) x (—1)^+^ at the sites immediately adjacent to the domain bound- 
ary, that is, at i x = 51 — i s , i y = 1 + i s with integer z s 's, in 8 steps from near zero to the full 
distortion close to s . At each step, we minimize the total energy with respect to the distortions at 
all other sites using the Euler method. This gives rise to the distortion profiles along the horizontal 



direction shown in Fig. 10 The 2D configurations for the red, green, and purple dots in Fig. 10 
are also shown in Figs. 11 ^a), [TT|b), and [TT|c), respectively, where the color represents s x . The 
results show that s x (i x , i y ) x (— at (51, 1) and (50, 1) grow together, compensating s x dis- 
tortions with opposite signs at the two neighboring sites, and the domain boundary advances by 
two interatomic distances. 

We define the effective location of the domain boundary, ddb, according to 

s*. - 0.040 

4 = 2x — x - , (59) 

0.308 - 0.040' V 



where s* represents the value of s x (i X} i y ) x (— chosen at i x = 51 in Fig. 10 and 0.040 and 



0.308 the values of s* before and after the domain wall moves by two interatomic distances. The 



minimized total energy E' tot is plotted in Fig. 12 for < < 2, which shows the energy barrier. 



We compare three energies, E\ = —2.946, E 2 = —1.668, and E 3 = —3.693 for ddb = 0.0, 1.0, and 
2.0, which correspond to the configurations shown in Figs. 11 a), [TT^b), and [TT|c), respectively. 



Most changes in distortion occur in the 64 x 2 plaquettes near the domain boundary, as shown in 



Figs. [10] and [TTJ Therefore, the energy difference between the two stable domain configurations 
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for d db = 0.0 and 2.0 is (E\ — E 3 ) / 128 = 0.0058 per site, which agrees with the energy difference 
per site, 0.0058, between the undistorted and distorted uniform phases in Fig. [3jb). The energy 
barrier normalized for 64 x 2 plaquettes, that is, (E 2 — £'i)/128, is 0.0100, which is of the same 
order of magnitude as the height of the energy barrier, AE = 0.0160, between the local and global 
energy minima in Fig.^b). From this analysis, the energy barrier against the uniform shift of the 
domain wall would be of the order of 2AE multiplied by the domain wall length in the units of 
interatomic distance, which would be a macroscopic energy barrier for the domain walls of micron 
length scale. We emphasize that discreteness of the lattice in our model is essential for this energy 
barrier, which is an example of the Peierls-Nabarro barrier.^ 



B. Stability against non-uniform domain wall modifications 

The importance of the long range interaction between strain fields is even more evident for the 
stability against nonuniform modification of domain walls. As an example, we convert a patch 
of the undistorted region in a configuration similar to Fig. [8ja), into a distorted state initially and 
then relax the whole lattice according to the Euler method. The initial, two intermediate, and fi- 



nal configurations are displayed in Fig. 13 The results show that the distortion in the converted 
region disappears initially except for two atomic layers [Fig. [T3fc)], which shrink laterally by 
further relaxation, restoring the original configuration [Fig. [T3~fcl)]. The simulation demonstrates 
the stability of domain structure against non-uniform modification of the domain walls. To gain 
further insight into the role of lattice compatibility, we examine other modes and energy distribu- 
tions. Figures [147a), [T4|b) and [T4fc) show the modes, ex, e 2 , and e 3 for the s x distortion given 
in Fig. |T3jb). The strain field tends to spread into the domains from the domain boundary. In 



particular, the field inside the converted patch in Fig. 14 c) cannot reach — C^sl/A^, the full 



distortion of e 3 inside the domain, due to the strain compatibility. The map of E to t(i), the sum 



of the terms with the site index i in Eqs. (39)-(41 ), is shown in Fig. 14 cl), which implies that the 



energy cost for creating the distorted patch is not confined immediately around the interface, but is 
distributed over the whole converted patch. This is different from systems with short range inter- 
actions only, for which the energy cost would be confined near the domain wall within the range 
of the interaction. This difference shows that the lattice constraint, leading to the effective long 
range interaction, plays an important role in the stability of phase coexistence against non-uniform 
modification of the domain boundary. Similarly, we find that if we convert a patch of the distorted 
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region of a similar size as above into the undistorted phase, the system relaxes back to the original 
configuration. 

However, the above results do not mean that it is impossible to convert a region between phases. 
For example, if we convert a large enough patch, as shown in Fig.[T5j^a), even though the distortion 
in the most converted region disappears initially, the two distorted layers remaining in Fig. \T5[b) 
expand laterally, as shown in Fig. [T5|c). Eventually, the distorted domain grows by two atomic 
layers, as shown in Fig.[T5|d). These results, particularly the different relaxation behavior for the 



configurations in Figs. \l3[c) and[T5[c), show that the energy barrier for the growth of the low 
energy phase involves simultaneous distortions of a significant number of unit cells just next to the 
domain wall. Slow growth of the low energy phase has been observed in a number of experiments 
for manganites. For example, Ref. Unreports a time scale of the order of 10 minutes for the growth 
of the low energy phase. A rough order of magnitude estimate of the energy barrier AE can be 
made by assuming an activated thermal process so that the relaxation time r = r exp(A_E//c B T), 
where t represents the intrinsic time scale for ion motion, ks is the Boltzmann constant, and T 
is temperature. With r ~ 10 3 s, r ~ 10~ 13 s, ksT ~ 10 meV, we obtain AE of the order of 
1 eV. If we consider the typical energy scale for the distortion of unit cell to be 1 - 10 meV, this 
energy barrier corresponds to about 100 - 1000 unit cell distortions within the layer just next to 
the 2-dimensional interface, consistent with our simulations. A similar growth of the undistorted 
region occurs if we convert a large enough distorted region into the undistorted phase, as shown 



in Fig. 16 The result in Fig. 16 is reminiscent of experiments in which the volume fraction of the 
undistorted phase is increased by external perturbations such as x-rays or lightpSH 



IV. SUMMARY 

We have discussed various aspects of a model for the strain-induced phase coexistence observed 
in perovskite manganites. A square lattice and associated atomic scale distortion modes were used 
to construct an energy expression with local and global energy minimum states, which captures 
features of manganites essential for phase coexistence: a local energy minimum metallic state 
without lattice distortions and a global energy minimum insulating state with short wavelength 
and uniform lattice distortions. Explicit expressions for modes, constraint equations, energies, 
and energy gradients have been presented. Our simulations for an energy landscape with a low 
energy barrier against transforming from undistorted local to distorted global energy minimum 
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states revealed nucleation with anisotropic correlation upon rapid quenching. Our simulations for 
an energy landscape with a high energy barrier showed stable coexistence of undistorted metallic 
and distorted insulating domains. Further, we studied the stability of such metal-insulator domain 
structures against various perturbations. We found that domain configurations are stable against 
uniform motion of the boundary due to the discreteness of the lattice and the intrinsic energy 
barrier between local and global energy minimum states. We expect that this intrinsic atomic scale 
energy barrier, multiplied by the number of atoms within the mesoscopic scale domain wall, is 
large enough to prevent the uniform motion of domain walls. For non-uniform modification of 
these walls, the long range interaction between strain fields gives rise to the domain wall energy 
distributed over the whole modified area for our 2D model (or volume for 3D system), rather than 
just the region confined near the domain wall, providing extra stability to the domain structure. To 
provide comparison, we carried out simulations for a system with a short range interaction only, 
which show no anisotropic nucleation or stable coexistence of local and global energy minimum 
phases. The above results demonstrate that the long range interaction between stain fields, and 
associated complex energy landscape, play an important role in metal-insulator coexistence in 
perovskite manganites. 

Establishment of more concrete connections between our model and experiments would be the 
goal of future studies. For example, the density of states at the Fermi energy level in the inho- 
mogeneous state can be compared with the conductivity measured in experiments. The effect of 
substrate-induced strain can be simulated in our model with additional energy terms representing 
the bonding between atoms in the film and the substrate. Thermal fluctuation can be simulated by 
the Monte Carlo method, which may provide insights into the origin of "strain" glass behavior that 
has been experimentally proposed as intrinsic rather than extrinsic.^ Furthermore, although our 
framework is based on the assumption that electronic and magnetic effects are adiabatically slaved 
to lattice distortions, our work can, in principle, be generalized to include these functionalities 
in a self-consistent manner. Such coupled models will be computationally intensive and our ap- 
proach has been to seek a minimal model. However, discrete strain or pseudo-spin models 43 with 
long-range interactions and disorder provide a skeletal approach to couple with magnetic spins and 
electronic densities within a mean-field or Monte Carlo scheme. Here the abundant literature on 
spin models is an advantage because even glassy behavior in electronic materials may be identified 
by an appropriate order-parameter. 
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Appendix A: Energy expressions for inhomogeneous states 

First, we represent ex, e 2 , and e 3 in the reciprocal space, and rewrite Ei and E c in Eqs. <|40|) and 



(41 ) in the following form: 



E i = N " E ^i(k)ex(-k) + ^e 2 (k)e 2 (-k) + ^e 3 (k)e 3 (-k), 



(Al) 



s y & 2 ]^(k)e^ 



(A2) 



Next, because the constraint equations apply differently depending on whether either k x or k y is 
zero or not, we divide the A;-sum into four parts, 

£=£ + £ + £ + £. (A3) 

% k x ^0,k y ^0 k x =0,k y ^0 h x ^0,k y =0 k x =0,k y =0 

and treat each of them separately. If k x ^ and k y ^ 0, constraint equations Eqs. (25)-p7|) are 
rewritten in the following way, which expresses the modes ex(k), e 2 (k), es(k) in terms of s x (k) 
and s y (k): 

i 

"71 



{ex)k x ^Q,k y ^o 
{e 2 )k x ^o,k y ^o 



k k - ' 

cot ~^s x (k) + cot Y s v( k ) 



(e 3 ) 



k x ^0,k y ^0 



71 

i 

71 L 



k k -> ' 

cot ~^s x {k) + cot -^s y {k) 

k k 
cot -^-s x (k) — cot -^-Sy(k) 



(A4) 
(A5) 
(A6) 



Therefore, the part with fc^O and k y ^ for E[ in Eq. ( Al ) is expressed as 



r 



(E[) kx ^0,ky^0 



AT2 



£ 

k x J=0,ky^0 \ s y 



B xx 


B xy 




' s x \ 


B xy 


Byy J 


k 


\s y ) 



where 



B xx (k) = \{M + A 3 ) cot 2 -j- + l -A 2 cot 2 ^ 
B yy (k) = \{Ax + A 3 ) cot 2 ^ + l -A 2 cot 2 ^, 
B xy (k) = l(A 1 + A 2 - A 3 ) cot y cot ^. 



(A7) 

(A8) 
(A9) 
(A10) 
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Similarly, the part with k x ^ and k y ^ for E c in Eq. (|A2b is equivalent to 



i k x ^0,k y ^0 



I 

71 



k k 
cot -^s x (k) - cot -^s y {k) 



J.k-% 



(All) 



For the terms with k x = and 7^ 0, we apply the constraint equation e\(k) + e 3 (k) = to 
eliminate e\(k) in E\ in Eq. (Al ) and obtain 

(Ei + E c ) kx=0 ^ = N 2 Y \( A ^ + A 3 )e 3 (-k)e 3 (k) + -A 2 e 2 (-k)e 2 (k) 

k x =0,k y ^0 

+ YC 3 [s x (z) 2 -s y (j:) 2 ] Yl e 3 (k)e^. (A12) 

Since we are interested in metastable phases in this work and e 2 (k) and e 3 (k) are independent for 
k x = and k y ^ 0, we minimize the energy (Ei 4~ E c ) kx= Q^ ky ^Q with respect to e 2 (k) and c 3 (k) 
separately and obtain 

C 3 



(ex) 
(e 2 ) 



mm 
k x =0,k y ^0 



k x =0,k y ^0 



_( \min 

0. 



Al + A; 



where 



?M-4) = ]!pJ2[ s ^) 2 - s y& 



-ik-i 



(A13) 
(A14) 

(A15) 



The minimized energy expression for (E t + E c ) kx=0:ky7 L is 

a 



(E t + Ec)™^^ = - I N 2 J2 ^-k y (s 2 x -s 2 y )^ ky (s 2 x -s 2 y ). (A16) 

Z{A X + A 3 ) 



k x =0,k y ^0 



We apply a similar analysis for the terms with k x 7^ and k y = in Eqs. (Al ) and (A2). Using 
the constraint e^k) — e 3 (k) = 0, we eliminate e\(k) and obtain 

(Ei + E c ) kx ^ ky=0 = N 2 \( A ^ + A 3 )e 3 (-k)e 3 (k) + -A 2 e 2 (-k)e 2 (k) 

k x ^0,k y =0 

+ ^C 3 [ S ,« 2 -^) 2 ] Y <^<' >r "- (A17) 



Separate minimization of this energy with respect to e 2 (k) and e3(/c) leads to 

C 3 



(ex) 



min 

kx ^~0)^T/' 



=0 = ( e s) 



^^0,^=0 _|_ ^ 



(^2)^0,^=0 — 0) 



(A18) 
(A19) 



(Ei + E c )™% >ky=0 = - °l N 2 Y -4)^.o(4-4)- (A20) 



k x =£0,k y =0 
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The terms with k = in Eqs. (All and ( A2 ) are 



(Ei + E c ) k 



x — Oj/Sy — 



N 2 



(A21) 



We minimize the above expression with respect to e±(k — 0), e 2 (k — 0), and e 3 (k = 0) indepen- 
dently, since they are not constrained to each other, and obtain 



( e i)™=0,fey=0 _ 

(e2)™= ,fc y =o = 

/ \min 
\ c 3)k x =0,k y =0 



-9lj. (s 2 -s 2 ) 



(Ei + E c 



2A- 



N 2 [^ k=0 (4 ~ * 



r 



(A22) 
(A23) 
(A24) 

(A25) 



Finally, by adding the terms with different cases of k x and k y found above, we obtain the following 
total energy, E' tot , which depends only on s x and s y : 



E Lt( s x, Sy) = E' l+C + E s 



(A26) 



where 



E' l+C = (Ei) kx ^o,k y ^o + (E c ) kx7 L ,k y ^o + (Ei + E c )™^ hy=0 + (Ei + E, 



mm 
cJk x =0,k y ^0 



{E l + E l 



(A27) 



and E s is given by Eq. (39). We use this energy expression E' tot (s x , s y ) for the simulations of 
inhomogeneous states. 

In addition to s x (i) and s y (i) configurations, e 1 (i), e 2 (i), and e 3 (i) configurations give useful 
information on the nature of the inhomogeneous states. The relations used to eliminate e l5 e 2 , 
and e 3 variables above, namely, Eqs. (|A4])-((A6]), ( |AT3| ), dAT4] ), ( |AT8] ), dAT9l ) and (|A22])-((A24]) 
for different cases of k x and k y , are used to find ei(k), e 2 (k), and e^(k) from given s x (i) and 



s y (i), which lead to ei(i), e 2 (i), and e 3 (i) configurations. Equations (28 )-(37 ) are used to find the 
displacements, u x (i) and u y (i), from the distortion modes. 
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Appendix B: Gradients of the energy expression for simulations using Euler method 

The gradient of E' tot (s X) s y ) necessary for the Euler method is found from 



dE' tot (s x ,s y ) = d(Ei) kxli0ikyji0 d(E c ) kx ^ 0iky ^ d{E l + E^ 1 ^^ 
ds x (i) ds x (i) ds x (i) ds x (i) 

| 8(E l + gjggp^ | d(E l + gjggp^ | d E s 



ds x {i) ds x {i) ds x (i) 



dE' tot (s x ,s y ) 

dSy(l) 



d{E{) Wy ^ + 9(£ c ) w + dim + E, 



mm 
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ds y (i) 

d(E t + E c )™o,k y =o , d(E l + E c j' k ^ ky=0 0E„ 



dSy(t) 



ds y (i) 

+ 



dSy(l) 



The expression for each term is given below. 



d{E t ) 



k x ^0,k y ^0 
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d(E c ) ^0,ky^O 

dSy(l) 
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FIG. 1. Two-dimensional square lattice with a monatomic basis, ^represents the displacement of the atom 
at the site with the index i = (i x , i y ), where integer i x and i y range from 1 to N. The site at the bottom 
left corner is chosen for i = (1,1). In this work, the lattice constant a is irrelevant for the expressions of 
modes and energies, and can be chosen arbitrarily, so long as the sequence of the atoms is not changed by 
the displacements. For figures in this work, a is chosen as 10, but the result itself is independent of the 
choice of a. 




FIG. 2. Symmetry distortion modes^for the motif for the two-dimensional square lattice with a monatomic 
basis shown in Fig.[T] 
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FIG. 3. (Color online) Energy landscape for the homogeneous phase, (a) Solid dots mark the locations 



of local and global energy minima for £^™ n /N 2 in Eq. (49) in the s x — s y plane, where s x = s x (tt, tt) 



and s y = s y (7T, tt). The local energy minimum states with the opposite signs of s x or s y are related by a 
phase difference of the short wavelength distortion, whereas the local energy minimum states along s x axis 
is related to those along s y axis by the exchange of x and y axes, (b) E^™™ /N 2 versus s x with s y = for 
the two chosen parameter sets. The parameter values for the deep local energy minimum case, represented 
by the thick red curve, are A x = 7, A 2 = 4, A 3 = 6, B = 2, C 3 = 20, G x = 60, G 2 = 80, Hi = 480, 
and H 2 = 640, for which G[ = -73.3, G' 2 = G 2 + 2C 2 /A 3 = 213, s = 0.34, and |e 3 | = 0.39, 
and E t Q™ m /N 2 has a value of -0.0058 for the global energy minimum states. The height of the energy 
barrier between local and global energy minima is 0.0160, measured from the local energy minimum. The 
only different parameter value for the shallow local minimum case, represented by the thin blue curve, is 
B = 0.5, which gives rise to so = 0.38, \e$\ = 0.49, the global energy minimum E\^ l%n /N 2 of -0.1053, 
and energy barrier of 0.0009. Such difference in the energy landscape can be related, for example, to the 
average size of the rare earth and alkali metal elements in manganites. 
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FIG. 4. (Color online) Simulation of lattice relaxation for a 32 x 32 lattice for the energy landscape with 
a shallow local minimum, namely the thin blue curve in Fig. |3jb). The vertices correspond to the positions 
of atoms and actual distortions are shown. The color represents p% = s% — Sy, green corresponding to 
zero, red and blue ±Sq, except ±(2.6so) 2 for the panel (a). The dynamics is governed by the Euler method, 
simulating a rapid quenching from a random initial configuration shown in (a), (b)-(g) show intermediate 
configurations and (h) the final stable configuration, in which most region is changed to distorted state 
except the antiphase boundaries represented in green. Highly anisotropic nucleation and the effect of long 
range interaction can be identified, for example, in (d) and (e). 
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FIG. 5. (Color online) Simulation of relaxation for a system with a short range interaction only. The 
variable p = M 2 is plotted on a 32 x 32 lattice with a periodic boundary condition. Red represents 
p = 28.6Mq for (a) and p = Mq for (b)-(h), and green p = 0. A random initial configuration is shown in 
(a), (b)-(g) show intermediate configurations, and (h) the final stable configuration. The result shows that 
nucleation for short range interaction is nearly isotropic, unlike the case in Fig. [4] 




FIG. 6. (Color online) (a) Map of local electron DOS at E = calculated for the distortion shown in 
Fig.[4]T). Blue, green, and red correspond to zero, 0.25, and 0.5 state per site per unit energy, respectively, 
(b) Local DOS calculated at i = (17, 10), the center of an undistorted metallic region, and at i = (7, 9), 
center of a distorted insulating region. It shows coexistence of metallic and insulating regions. 
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FIG. 7. (Color online) Simulation of relaxation of lattice distortions for random initial distortions in a 
narrower range compared to the case in Fig.[4j Color scheme is identical to Fig. |4] The initial configuration 
in (a) is obtained by multiplying 0.5 to the random initial s x (i) and s y (i) in Fig.^a). The comparison with 
Fig. ^indicates that the low temperature phase depends on how the system is prepared, as observed in some 
manganites. 
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FIG. 8. (Color online) (a) Stable configuration of distorted and undistorted domains for a 64 x 64 lattice 
for the energy landscape with a deep local minimum, shown by the thick red curve in Fig. |3jb). The color 
represents j»3 with red and green for s^, and 0, respectively, (b) Map of the local electron DOS calculated at 
E = for the distortions in (a). Red, green, and blue correspond to 0.6, 0.3, and state per site per energy, 
respectively. There is no intrinsic length scale for the size of the domains in our model, and the size of the 
domain depends only on the predesigned initial condition and can be as large as micron in principle, which 
is consistent with experimental observations. 
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FIG. 9. (Color online) Local electron DOS calculated at the center of an undistorted region (thick green 
line), i = (36, 1) in Fig. [8ja), and at the center of a distorted region (thin red line), i = (6, 1) in Fig.[8[a). 
Unlike the case of small domains in Fig. [6] a clear gap opens inside the large domain of distorted phase. 
The local density of states is always symmetric with respect to E = 0. Therefore, for Ep = 0, the electron 
number is 0.5 at all sites, demonstrating that the metal-insulator domain structure does not requires charge 
inhomogeneity. 
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FIG. 10. (Color online) The red dots represent the profiles of s x (i x ,i y ) x (— with i y = 1 near 

the domain boundary in Fig.[8ja). Other dots show how this profile changes as the domain boundary shifts 
uniformly by two interatomic distances. Lines are drawn to guide eyes. The final configuration shown in 
purple dots is equivalent to a parallel shift of the initial configuration shown in red dots, but the intermediate 
configurations are not, due to the discreteness of the lattice. 
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(c) 




FIG. 1 1 . (Color online) Configuration of lattice distortions near the domain boundary, as the boundary 
moves by two interatomic distances. The color scheme is different from Figs.|4j|7j and[8ja). Here, the color 
represents the s x mode, with red, green, and blue corresponding to s a , 0, and — s a , respectively, (a), (b), and 



(c) correspond to the profile represented by red, green, and purple dots in Fig. 10 Configuration (c) has a 
lower energy than (a), but the higher energy of the intermediate configuration (b) serves as a Peierls-Nabarro 



barrier, as shown in Fig. 12 
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FIG. 12. Total energy E' tot for the 64 x 64 system, given by Eq. ( A26 1, versus the location of the domain 



boundary defined by Eq. (59 1, as the boundary moves by two interatomic distances. Each point is found 
from each corresponding curve in Fig.[l0j The configurations in Figs.[TTJa),[TTJb), and[TTJc) correspond to 
ddb = 0.0, 1.0, and 2.0, respectively. The energy barrier prevents a parallel shift of the domain wall. 
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FIG. 13. (Color online) Simulation of the domain wall stability against a small nonuniform modification of 



the domain boundary for the configuration similar to Fig. |8[a). The color scheme is identical to Fig. 1 1 (a) 
represents the initial perturbed configuration, (b) and (c) show intermediate configurations, (d) represents 
the final stable configuration, which is identical to the original configuration before the perturbation. It 
shows that the original domain configuration, such as the one shown in Fig. [8] is stable against small non- 
parallel shift of the domain boundary. 
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FIG. 14. (Color online) Colors in (a), (b) and (c) show the e\{i), e2{i), and es(i) for the configuration 



shown in Fig. 13 b). Colors in (d) show E to t(i), which is the sum of the terms with the site index i in 



Eqs. d39j)-(41). In the panels (a), (b), and (c), red and blue correspond to ±0.45 and green to zero. In 



the panel (d), red and blue correspond to 0.06 and -0.006, respectively. Typical values of e$(i) and E to t{i) 
inside the converted patch are -0.08 and 0.02, respectively. Panel (d) shows that the energy cost for changing 
domain wall in a non-parallel way spreads over the whole changed area (volume in 3D), which is different 
from systems with a short range interaction only and is responsible for the stability against such non-parallel 
shift of the domain walls. 
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FIG. 15. (Color online) Simulation of the domain wall shift by large enough nonuniform modification 
of the domain boundary, converting a large undistorted patch into a distorted state, for the configuration 
similar to Fig. [8f a). The color scheme is identical to Fig. [TT] An initial configuration is shown in (a), (b) 
and (c) show intermediate configurations, (d) represents the final stable configuration, which shows that the 
distorted region has expanded by two atomic layers. It shows that if the perturbation is large enough, the 



system relaxes to a different metastable configuration, unlike the case in Fig. 13 
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FIG. 16. (Color online) Simulation of the domain wall shift by a large enough nonuniform modification 
of the domain boundary, converting a large distorted patch into an undistorted state, for the configuration 
similar to Fig. [8ja). The color scheme is identical to Fig. [TT] An initial configuration is shown in (a), (b) 
and (c) show intermediate configurations, (d) represents the final stable configuration, which shows that 
the undistorted metallic region has expanded. This simulation may be related to the experiment results in 
Ref. HH in which the x-rays destroy the charge-orbital ordering and the short wavelength distortions and 
increase conductivity. 
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